Analyzing topographic change in the USA

Check out the data

topo_change <- st_read('data/topo_change/topo_change_polygons.shp')
## Reading layer `topo_change_polygons' from data source `C:\Users\jesca\OneDrive\Documents\WR674\mining_impacts_to_elevation\data\topo_change\topo_change_polygons.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8958 features and 51 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -2332814 ymin: 407783 xmax: 2238422 ymax: 3128340
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#Checkout the column names
#names(topo_change)
biggest_area <- topo_change %>%
  arrange(desc(AREA_SQ_KM)) %>%
  slice(1:10)
# Map all mines
mapview(biggest_area)

Checkout mine in Arizona

# Subset to just Arizona
az_mine <- biggest_area %>%
  filter(QUADNAME == 'esperanza_mill_AZ')
#Check that it is the right site
mapview(az_mine)

Getting our own elevation datasets

Download pre-mining DEM

#Check projection of az_mine
#st_crs(az_mine)
az_raster_before <- get_elev_raster(az_mine,z=12)
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#Look at the structure of the data
#str(az_raster)
#Summary of the data
#summary(az_raster)
plot(az_raster_before)

#mapview(az_raster_before)

Shrink the resolution

small_az <- aggregate(az_raster_before,fact=3)
mapview(small_az)

Take the slope of the data

az_slope <- terrain(small_az,opt='slope',unit='degrees')
mapview(az_slope)

Take the aspect of the data

az_aspect <- terrain(small_az,opt='aspect',unit='degrees')
mapview(az_aspect)

3D plotting mechanism for elevation data

az_matrix <- matrix(raster::extract(small_az,
                    raster::extent(small_az),buffer=500),
                    nrow=ncol(small_az), ncol=nrow(small_az))
az_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(az_matrix),color='desert') %>%
  add_shadow(ray_shade(az_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(az_matrix)) %>%
  plot_3d(az_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)

In Class work

For this section we are really going to focus on the mine area so first I’m going to buffer our az_mine polygon by 2km and then clip both the before and after mining rasters to this new shape

#Read in raster data
az_raster_after <- raster('data/srtm_14_06.tif')
az_mine_2km <- az_mine %>%
  st_transform(2163) %>%
  st_buffer(2000) %>%
  #Reproject into SRTM datum (WGS84)
  st_transform(st_crs(az_raster_after))
#Get extent of buffered mine area
aoi <- extent(az_mine_2km)
#Crop after to area of interest
mine_after <- crop(az_raster_after,aoi)
#Crop before to area of interest (need to also reproject)
mine_before <- crop(az_raster_before %>%
                      projectRaster(.,crs=projection(mine_after)),
                    aoi) %>%
  projectRaster(.,mine_after) #match resolution
#Force rasters to be exactly same size
mine_before_same <- trim(mine_before)
mine_after_same <- crop(mine_after,mine_before_same)
#Stack these things in a rasterBrick (which we can do because we 
#made them the exact same size)
mines <- brick(mine_before_same,mine_after_same)
#Add a dem_difference layer
mines[[3]] <- mines[[1]] - mines[[2]]
#Rename raster layers
names(mines) <- c('Pre_mining','Post_mining','DoD')
plot(mines)

All data for this section should be the data clipped to the mining area (generated in the code chunk above)

1) Get the slope of the SRTM (post-mining) dataset and plot this data (use the data that is clipped to the mine only).

after_slope <- terrain(mine_after_same,opt='slope',unit='degrees')
plot(after_slope)

2) Use the brick code from above to plot slopes for the pre-mining, post-mining, and difference of slope

before_slope <- terrain(mines[[1]],opt = 'slope', unit = 'degrees')
plot(before_slope)

diff_slope <- terrain(mines[[3]], opt = 'slope', unit = 'degrees')
#plot(diff_slope)

mine_slopes <-brick(before_slope, after_slope, diff_slope)
names(mine_slopes) <- c('Pre_mining','Post_mining','DoD')
plot(mine_slopes)

3) Use rayshader to make a 3d plot of the pre-mining mine

pre_matrix <- matrix(raster::extract(mine_before_same,
                    raster::extent(mine_before_same),buffer=500),
                    nrow=ncol(mine_before_same), ncol=nrow(mine_before_same))
pre_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(pre_matrix),color='desert') %>%
  add_shadow(ray_shade(pre_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(pre_matrix)) %>%
  plot_3d(pre_matrix,zscale=10,fov = 0,theta=135,zoom=0.85)
render_snapshot(clear = TRUE)

4) Use rayshader to make 3d plot of the post-mining mine

post_matrix <- matrix(raster::extract(mine_after_same,
                    raster::extent(mine_after_same),buffer=500),
                    nrow=ncol(mine_after_same), ncol=nrow(mine_after_same))

dim(post_matrix)
## [1] 88 80
dim(pre_matrix)
## [1] 88 80
post_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(post_matrix),color='desert') %>%
  add_shadow(ray_shade(post_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(post_matrix)) %>%
  plot_3d(post_matrix,zscale=10,fov = 0,theta=135,zoom=0.85)
render_snapshot(clear = TRUE)